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Abstract 



A brief review of the Hadron Gas model with reference to high energy heavy 
ion collisions is presented. The entropy dependence on temperature and baryonic 
chemical potential is numerically calculated, together with the entropy distribution 
between baryons and mesons. The theoretical entropy for a QGP with equivalent 
parameters is also calculated. It is shown that at low temperatures the dominant 
entropy contribution comes from baryons, while at high energies the dominant 
contribution comes from mesons. The turnover from baryon to meson domination 
occurs at T 140MeV, which corresponds to an energy domain of the AGS at 
Brookhaven National Laboratory. 
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Chapter 1 
Introduction 



Relativistic heavy ion collisions provide us with the opportunity to study strongly 
interacting matter at extremely high temperature and pressure. It is the current 
goal of collision experiments to attempt to create a new state of matter known 
as a Quark-Gluon Plasma (QGP). This refers to a state where the quarks and 
gluons, which were originally bound into nucleons, become deconfined and form a 
self-consistent plasma. This form of matter is believed to have existed in the very 
early universe (shortly after the Big Bang). 

The Quark-Gluon Plasma was originally proposed as a byproduct of QCD. Quan- 
tum chromodynamics was developed as a theory of the strong interaction, but 
unlike QED, it is almost impossible to obtain analytic solutions. Lattice QCD 
(discrete points represent space-time) has been worked on numerically, and it sug- 
gests that a QGP forms when the temperature (~ 150 — 200 MeV) or density (up 
to 5-lOpoj Po ~ 3 x 10 14 g/cm 3 , the nuclear density) becomes large enough. 

The search for the QGP is currently ongoing. As the beam energies of heavy ion 
accelerators continue to increase, the systems resultant from the collisions have 
initial temperatures and densities close to the above thresholds. 

Since direct observation of the QGP is impossible (it is very short-lived, and its 
component quarks and gluons become confined into hadrons before they reach 
the detectors), one has to infer its formation from the observed distributions of 
hadrons, leptons and photons. This is highly non-trivial, as most hadronic ob- 
servables reflect conditions at freeze-out (when the hadrons stopped interacting - 
similar to the surface of last scattering.), by which time final state interactions 
may have caused loss of all information about the initial partonic state. 

After the release of the latest data from RHIC, there has been some speculation 
at to whether data showed signs of existence of QGP, but the issue has not been 
settled yet. It is hoped that with still higher energies expected in the near future 
at LHC, results signalling the formation of QGP will be observed. 
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1.1 Ideal Hadron Gas Model 



While QGP is expected to be the high-energy limit of a collection of hadronic 
particles, at lower temperatures there is a different model. The conventional, 
confined phase is usually referred to as Hadronic Gas (HG), consisting of hadrons 
of different types (including the short-lived resonances) such as 7r, p, N, A, K, 
whose properties (mass, spin, degeneracy) are well known. 

This phase does exist in heavy ion collisions. Even if a QGP is created, it will 
rapidly expand until the temperature drops below the critical value, and hadroniza- 
tion takes place. This means that the deconfined quarks and gluons of the QGP 
will form bound states of baryons, mesons and anti-baryons. 

The system will then form a small volume that is filled with a huge variety of 
interacting hadrons. Assuming that a state of equilibrium is reached it seems 
plausible to treat the system (fireball) using the methods of statistical mechanics. 
The Thermal Model can be applied to a system at thermal equilibrium, although 
to find the particle multiplicities, chemical equilibrium also required. 

1.1.1 Thermal and chemical freeze-out 

Directly after the hadronization, the particles will be sufficiently close to interact 
strongly (distances of order lfm), and it is expected that this form of interaction 
will cause the system to reach a state of "chemical equilibrium" (chemical referring 
to the composition of the fireball). Hence, chemical equilibrium means that the 
number of each form of particle will be constant - the rate of creation/annihilation 
for each type of particle will be exactly equal. This will then continue until the gas 
cools further (it is expanding), and the temperature drops below a critical value. 
This is called chemical freeze- out - the point at which inelastic collisions cease and 
all particle ratios are frozen. 

A second form of freeze-out, known as thermal freeze- out is often also introduced. 
After the chemical freeze-out, once particle ratios become fixed, particles still in- 
teract elastically (eg. 7r + A^A^7r + A). These types of reactions have larger 
cross-sections and continue to occur, redistributing momentum in the system while 
unchanging the chemical composition. The fireball continues expanding and cool- 
ing until these interactions cease (thermal freeze-out). The particles then fly off 
towards the detectors. 

By observing the particle abundances (which show the properties of the system 
at chemical freeze-out), one finds that (at CERN SPS energies) the distribution 
is characteristic of one at T ~ 170 MeV. By observing the particle momentum 
distributions (which characterize the system at thermal freeze-out), one can infer 
that (for CERN SPS), the thermal freeze-out occurs later at T ~ 130 MeV. 
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1.2 Theoretical Formulation of the Thermal Model 



In the fireball, the number of particles is not fixed, and so one works with a Grand 
Canonical Ensemble. Firstly, consider the analysis of a static fireball. Hence the 
system lives for a period of time in a fixed volume V. Then for one type of particle: 

Vt GC = -kT In Z (1.1) 

where Z is the partition function, the form of which depends on whether the 
particle follows Bose-Einstein or Fermi-Dirac statistics. The specific quantities 
will be discussed in detail in the next section. The multiplicity of hadrons of 
species % is given as: 



N t = Vru 



V 9i 



(1.2) 



where Hi is the number density, gi = 2Jj + 1 is the spin degeneracy factor, and 
fi(p) is the momentum distribution function. In a thermodynamic equilibrium 
distribution, these functions have a relatively simple form: 



Up) = 



kT 



+ e = 



yV 



+ e 



-i 



(1.3) 



Here k is Boltzmann's constant, T is the temperature of the system, and /ij is 
the chemical potential of hadron i. The quantity e equals +1 for fermion (FD 
statistics), and —1 for bosons (BE statistics). The limit e — > corresponds to 
classical (Boltzmann) statistics. 

The chemical potential Hi is what ultimately distinguishes the hadrons from each 
other. The quantities which are conserved for particle interactions within a hadron 
gas (which is sufficiently short-lived and short-range that it only interacts strongly) 
are the baryon number, charge and strangeness. The chemical potential is then a 
linear combination of the three potentials: 



(1.4) 



Here, Bj, Qj and S« are the baryon number, charge and strangeness of the i'th 
hadron respectively. This can be extended by adding the charm and bottomness 
potentials in a similar linear fashion. However, for the current calculations, only 
the particles composed of the u, d and s quarks have been used. 

Introduction of the chemical potentials hq and /is allows us to fulfill the 
appropriate conservation laws (for strong interactions). The strangeness of the 
system must be zero: 

J2^ = VJ2^ = (1.5) 



This can be implemented without knowing anything about the volume V. 
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Next, the total charge of the fireball must be the same as the total charge of the 
colliding nuclei. This calculation is however complex (volume modulated), but one 
needs to observe that the same factor exists in total baryon number. Dividing, we 
require the electric charge to baryon ratio to be: 



J2i B i N iJ HG \T,i B i N i J nuclei A 

This can be implemented for a general collision of two nuclei by varying \iq until 
the above equation is satisfied. For Pb+Pb and Au+Au collisions, the ratios differ 
very slightly, but can be very well approximated by: 

N Z 

» 0.66 and - w 0.40 (1.7) 

reflecting that there are slightly more neutrons than protons. A simplification I 
was encouraged to make was to to fix the potential fj,g — 0. As m, p ~ m n to very 
good accuracy, setting the electric charge potential to zero introduces a charge- 
independence into the system, and as protons and neutrons have almost the same 
mass, the model will produce an equal number of them. Hence 4 = 0.5 for \iq = 0. 

This may seem like an unnecessary simplification, and a deviation from real heavy 
ion collisions. However this is still exact for O+O and S+S collisions and it does 
produce a rather large simplification in the working. As it is, provided with T and 
fiB, there is only one unknown parameter (fis), and although there is no analytic 
method to solve for fi$ to satisfy (1.5), one can do it iteratively. Without specifying 
fiQ, there are two unknown parameters, and to satisfy both strangeness and charge 
conservation requires iterative methods in 2-D. 

If one does not implement this simplification and for a given T and [is does actually 
work out the value of \xq for a heavy ion collision, one finds that it is usually very 
small and negative (/iq < 0, to make neutrons more favourable than protons). 
Typical values are —10 MeV < \iq < — 4 MeV, depending on the temperature. 

Furthermore, the actual number of particles iVj and entropy Si depend on the 
volume V of the fireball, and for this project that is an unknown quantity. I 
therefore always worked with the number and entropy densities and Sj. However 
all the particle proportionalities are still known. 

The above is true for a static fireball - one that occupies a fixed volume V. In 
a more general case, when the expansion of the system at freeze-out cannot be 
neglected, one needs to use the more complex Cooper-Frye formula to calculate 
the total yield of particles. However it has been shown pQ that in this more complex 
system, the particle ratios stay the same as in the static fireball lone; as 

the thermodynamic parameters are constant along the freeze-out surface, so to get 
the particle ratios, one can work with a static fireball. 
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Chapter 2 

Numerical Integration 



The hadrons come in two varieties: baryons (anti-baryons) which obey Fermi- 
Dirac statistics, and mesons which obey Bose-Einstein statistics. Although the two 
behaviours are very different at low temperatures, the formalism is very similar - 
often just a change of sign. I'll thus deal with both cases using ± and =F notation, 
with convention (top sign for fermions). 



2.1 Quantum Statistics Integrals 



2.1.1 Number densities 



As mentioned in the previous section, we define a momentum distribution function: 

-l 



fi = fi(p) = 



e kT ±1 



e *r ± i 



Then the number density of hadron % is given by: 

d 3 p 



rii = gi 



fi(p) 



(2.1) 



(2.2) 



(27Tft) 3 ' 

Assuming that the fireball is isotropic, the integral separates into angular and 
spatial parts: 

"°° p 2 dp 



rii = 9, 



/2tv pi 
dxj> J d(cos9) 



r/i(p) 



9i 



fi(p)p 2 dp (2.3) 



/o (2nh) 3J ^ J 27r 2 n 3 J(] 
Going over to a system of units where H — 1, the multiplicity of hadron i becomes: 



rii = 



fi(p)p 2 dp 



9i 
2tt 2 



f 

Jo 



p 2 dp 



e kT ± 1 



(2.4) 
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2.1.2 Entropy densities 



For a set Grand Canonical ensemble of particles % in volume V, the partition 
function obeys: 

fip Wl±e-^ (2.5) 



lnZ = ± ^-/(^ k ( 1±e "^) 



Then the entropy is defined by: 

, , „ kTdZ 
= k\nZ + — — 

z or 



In lie kT 



{2txK) 



)) 



d 3 p 



111 lie kT 



= ±kVg 

= ±kVg. 

= ±kVg 

= ±kVg. 

= ±kVg 

= ±kVg 

= ±kVg 

= ±kVg 

= ±kVg 

= kVg t 



{2nhf 
d 3 p 

d 3 p 
(2nh) 3 

d 3 p 
{2-Khf 

d 3 p 
{2nhf 

d 3 p 
(2nh) 3 

d 3 p 
(2-kK) 3 

d 3 p 
(2nh) 3 

d 3 p 
{2-Khf 
d 3 p 



d 3 p d 
{2nh) 3 &f 

ln(l±e - ^) + / 



lie kT 



In 



± 1 

-JH 



± 



e kT 



E - Hi 
kT 



e fey ±1 



In I e ^ \ -\-\n\e ^ ± 1 j ± 



E - Hi 



kT 

-ln/i + ln 
-ln/i + ln 



E- fi. 



E-Hi 



kT 



kT 
-1 ± /i) 



g — A*i 
kT 



fi 



e fcT ±1 



e 



e fcr ±1 



(1 =F fi) 



(e^±l) (-l±/0 



:-i±/i) 



[- In /, - ln(l =F /,) (1 =f /,) + ln(/,) (1 T fi)} 
frfilnfi-Hl^mi^fi)} 



-fi\nfi T (lTfi)HlTfi)} 



(2nh) 3 

Integrating through the angular parts, one obtains the entropy density (h — k — 1): 



Si = 



2tt 2 



POO 

/ P 2 dp[-fikifiT(lTfi)HlTfi)} (2-6) 
./o 



2.1.3 Change of variables 



The integrals for n, and Sj are in p, but for p^> m: 

-l 



kip) 



± 1 



e t ± 1 



-1 






e T 



e r 



(2.7) 



So the integrand is modulated by a decreasing exponential of form e~ p ^ T . This 
is not desirable, as one would like the integrals to not be dominated by terms 
involving the temperature (to be able to observe qualitative behaviour). This 
is easily fixed by introducing x = p/T, =>■ p 2 dp = T 3 x 2 dx, and the limits of 
integration remain to oo. The two densities are then: 



n. 



2tt 2 

T 3_9i_ 

where we now have: 



j.3 9i i i: . ,., ,2 



fi(x)x dx (2. 
[-/i(a;) In /((z) =F (1 =F fi(x)) ln(l =f x 2 rfx 



1 



(2.9) 



The quantities m^/T and /Xi/T are constants for each particle and parameter set. 



2.2 Gauss-Laguerre Integration Technique 

The above integrals can be evaluated using a standard Simpson method of dividing 
the domain into strips and approximating the function on each strip by a parabola; 
the total integral being the sum of the areas of all the strips. This method is very 
accurate (provided the strip widths are small), but also time-consuming. I did 
implement this method, integrating the domain from to 20, with Ax = 0.002. 
However there exists a more efficient and far more subtle integration technique. 



2.2.1 Gaussian quadratures 

The method of Gaussian quadratures is not a very well known method of calcu- 
lating integrals, as it is only readily applicable to a limited set of integrands. This 
topic can be found in the literature, although the proofs of relevant theorems are 
not presented, and I've been unable to reproduce them. The following arguments 
are taken from Numerical Recipes combined with my own working. 

Consider the class of functions that take the form of a polynomial multiplied by 
some known function W(x). Then given the function W(x) and some integer N, 



9 



one can create a set of iV weights Wi and abscissas Xi, to approximate the integral 
by a finite summation: 

r b N 

/ W(x)f(x)dx^^2wif(xi) (2.10) 
Ja i=i 

such that the approximation is exact if f(x) is a polynomial. Now for a given N, 
there are 2N parameters that can be chosen at will (the weights and abscissas). It 
then follows that ()2.10|) can only be exact for 2N linearly independent polynomials, 
and it is usually chosen such that it holds for all f(x) with degree D < 2N — 1. 
For such a choice, approximating a more general integral by a finite sum is only 
a accurate if f(x) can be "well approximated by a polynomial". However almost 
all well-behaved (infinitely-differentiable) functions can be well approximated by 
a polynomial of sufficiently high degree. The value of N ought to be chosen 
sufficiently large to accommodate for slightly more complicated functions. 



Orthogonal polynomials 

For a specified function W(x) and limits of integration, define a "scalar product": 

(f\g) = [ W(x)f{x)g(x)dx (2.11) 

J a 

Then one can find a set of polynomials that satisfy (i) there exists exactly one 
polynomial of order j, called (pj(x)), and (ii) all the polynomials are mutually 
orthogonal over the specified weight function W(x). As a further requirement one 
can require that the polynomials Pj(x) be normalized - the scalar product with 
themselves giving unity. Then: 

(pi(x)\ Pj (x)) = 6a (2.12) 

A set of such orthonormal polynomials can constructed through a recurrence re- 
lation (Grand-Schmidt procedure), although it may not be the most efficient way. 

It can be shown that the polynomial Pj(x) has exactly j distinct roots in the 
interval (a,b). This becomes relevant when one uses the fundamental theorem 
of Gaussian quadratures: The abscissas Xi of the iV-point Gaussian quadrature 
formula ()2.10j) with weighting function W(x) on (a, b) are precisely the roots of 
the orthogonal polynomial p^{x). 

Once the abscissas known, the weights Wi can be found by solving: 

N „b rb 

^2p k (xi)wi = / W(x)p k (x)dx = / W(x)p k (x)dx ■ 8 0k (2.13) 

i=l J a J a 

for each value of k from to N — 1 (all integrands except for k = are zero, due 
to orthogonality). The parameters x^ and will make ()2.10|) exact for all any 
function that can be written as a linear combination of Pi(x). 
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2.2.2 Laguerre polynomials 



The form of integrals that are useful for this problem are polynomials modulated 
by e~~ x , integrated over the positive real axis. For this purpose take W(x) = e~ x , 
and the limits of integration from to oo. The orthogonal polynomials can be 
defined in terms of a recurrence relation, although there is a simpler method: 

• Define L (x) = 1. Then J °° L (x)W(x) = 1 as required. 

• For all other n define L n = ^-j-A^(x n e~ x ), from which one gets: 



d 



L n (x)W(x)dx = j -.^-{x n e- x )dx = - 

n\ dx n n\ 



(x n e- x ) 



dx 71 ' 1 Jn 



Directly showing that ()2.1Hj) is satisfied, and also giving the first orthonormality 
condition (as L = 1), if one rewrites the equation as J °° L (x)L n (x)W(x) = 5 0n . 

To see the form of the Laguerre polynomials, the first four have been calculated: 

L (x) = 1 
L>i(x) = 1 — x 



L 2 (x) = \(2-Ax + x 2 ) 
L 3 (x) = i(6-18x + 9x 2 



x 3 ) 



A general formula can be guessed, and shown to be true. To find the coefficient 
of x k from the expression -^(x n e~ x ), the product rule is used repetitively, and 
the derivative must be applied k times to the exponential and n — k times to the 
monomial. This produces a factor of (— l) k ■ |f. As there are n derivatives applied 
in total, and k are to the exponential, there is an additional factor of Q). Hence: 

Ln{x) 4f(-l)^f"V" = yVl) fe 7 %7T7^x k (2-14) 

v ' n\^ y ' k\\k ^ y ' (n - k)\(k\) 2 V ' 



Orthonormality 

The orthonormality of Laguerre polynomials is mentioned in all the literature, 
although I have been unable to find a proper proof of this property. The nor- 
malization J °° L\ (x)W (x)dx = 1 can be derived from (|2.14j) by calculating the 
coefficient of each x k , and using the property J °° x k e~ x dx = kl. This method is 
however very long and algebraically intensive, and not shown in this project. To 
show that distinct Laguerre polynomials are orthogonal is a harder task, one that 
I've been unable to complete myself, or indeed find a proof of this relation. 

If one takes it for granted that Laguerre polynomials are orthonormal, then one 
can expand any polynomial f(x) of degree less that 2N as a linear combination of 
Li(x), and as integrals preserve associativity, ()2.1()|) will be exact for all such f(x). 
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2.2.3 An explicit construction 

For iV = 2 it is possible to construct the abscissas and weights directly, by simply 
assuming (j2.1()j) and J °° x k e~ x dx = k\. This gives a system of four equations: 

Wi + w 2 = 1 

+ X 2 W 2 = 1 

x\wi + x\w 2 = 2 
xfwi + xlw 2 = 6 

Multiplying the second equation by {x\ + x 2 ) yields: 

Xi+ X 2 = (Xi + X 2 )(x 1 W 1 + X 2 W 2 ) = x\wi + XiX 2 (Wi + W 2 ) + x\w 2 = 2 + XiX 2 

While multiplying the third equation by (x\ + x 2 ) yields: 

2(x 1 +X 2 ) = (Xi + X 2 )(x\wi + X 2 2 W 2 ) = x\wi + XiX 2 (XiWi+ X 2 W 2 ) + x\w 2 = Q + XxX 2 

giving two equations for (x\+x 2 ) and Xix 2 . It follows that x±x 2 = 2 and x±+x 2 = 4. 
Solving for the individual x's (ordering is not defined - take x± < x 2 ) gives 

£i = 2 - V2 and x 2 = 2 + V2 

Note that the x^s are roots of the quadratic x 2 — Ax + 2, which is (up to a constant 
factor) L 2 (x). Hence the abscissae have been shown to correspond to the zeroes 
of the relevant Laguerre polynomial. 

The problem is now essentially solved. Putting the found values for the x^s into 
the original system leaves a set of linear equations for w^s, which can be trivially 
solved. The solution: 

w 1 = -(2 + V2) and w 2 = i(2 - y/2) 

A little algebra will show that the integral equations for (j2.10j) are satisfied for x k 
for k = 0,1,2,3. Hence the two Laguerre point problem is solved. This result, 
although simple to derive algebraically, allows one to integrate (from zero to in- 
finity) any exponentially modulated polynomial up to cubic by simply evaluating 
the function at two points. 

A result not really evident for the above case is the following. Suppose we use N 
Laguerre points. Then ANY polynomial P(x) which has zeroes at the Laguerre 
points will have a zero integral. This despite the extra N — 1 supposed degrees of 
freedom. An extreme case is a polynomial of the form: 

P{x) = {x N - x) x n^ x (x - x t ) 2 

which is positive for all x < xn (without loss of generality xn > Xi), and only 
negative once x > x^, by which time P(x)e~ x is expected to be heavily suppressed. 
This is not actually true: P(x) grows as x 2N_1 at x ~ x^, which is much faster 
than the negative exponential at the same point (for iV = 15, x^ ~ 48 3> e XN ^ 2N ^). 
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2.2.4 Applicability of Laguerre integration 



Unfortunately, this neat form of integration is only exactly applicable to integrals 
of the form (2.10) - polynomials modulated by a decreasing exponential. However 
the integrals that are needed for quantum statistics, are not like that (2.9). For 
large values of x, fi(x) — > e~ x . However, there is the e = ±1 factor, which creates 
slight problems for x ~ 0, and also the TOj/T factor which deforms the exponential 
factor. 

The method can however be salvaged by taking a "large" number of points. For 
n — 15 (as used), the method works for any polynomial up to degree 29. Although 
I have no analytic proof, it is rather intuitive that the method should still give a 
good approximate answer. There are 15 points where the function is evaluated. 
The integral is completely determined by the 15 data points. Now consider an 
exponentially modulated polynomial function passing through the same 15 points. 
There are still 14 degrees of polynomial freedom, by using which one can approx- 
imate fi(x) very well over it's domain. 

Numerical results for performing the integrals using the Gauss-Laguerre and Simp- 
son method (considered almost exact) show a difference between the two methods 
of less than 5 ppm (parts per million). For more complex functions (such as Sj), 
the Laguerre method appears even more doomed, as Sj involves terms proportional 
to (1 =f fi) ln(l =f fi), where the exponential fall-off is not clear. 

Consider the part of the expression for a; 3> ^. Then the function / of form 

(2.9) can be reasonably well approximated by: fi(x) ~ \e x ~^/ T ± l] 1 <C 1. Then 
using the Taylor expansion for the natural logarithm ln(l + x) ~ x, the entropy 
integrand becomes: 

-/iln/ iT (l=F/i)Ml=F/i) » -/iWi=F(l=F/0(=F/i) 

= fi(i^fi -hi ft) 

~ /i(l-ln/0 

« (1 + (x - fi/T))e ll/T - x 

= e" /T (x + l- n/T)e~ x (2.15) 

where I have approximated the distribution function by the Maxwellian limit for 
large x. So the entropy integrand, although complicated for x <C m/T, is in 
the required form of a polynomial modulated by an exponential for large x. So 
Laguerre integration can be applied here as well, although more cautiously than 
before. The integrals remain accurate (for the worst data set) to within a half of 
a percent. 
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Chapter 3 

Entropy Calculations 



The derived equation ()2.6|) gives an exact value for the entropy density of hadron i. 
To calculate the total entropy density of the system, one is required to sum over 
the contributions of all the hadrons. To this end I was provided with a list of all 
known hadrons composed of the u, d and s quarks that have a mass below 2.6GeV. 
This contains all the unflavoured and strange mesons, baryons and antibaryons. 
Charmed particles were not included, but the program can be easily extended 
to cater for these additional hadrons. The particle listing was presented in the 
following format (corresponding to 7T + , p°, K°, p, A ++ ) 

1. 0.140-1. 0. 0. 1. 

3. 0.770 -1. 0. 0. 0. 

1. 0.498 -1. 1. 0. 0. 

2. 0.938 1. 0. 1. 1. 

4. 1.232 1. 0. 1. 2. 

The first column corresponds to the degeneracy factor g s — 2s + 1. The second 
column gives the mass of the hadron in GeV. The third column gives e = ±1, 
depending on whether the particle obeys Fermi-Dirac or Bose-Einstein statistics. 
The next three columns give the strangeness S, the baryon number B and the 
charge Q of the particle. These three parameters will be used in the calculation 
of the particle's chemical potential (jl.4j) . These are all the parameters that are 
required to be able to work out the entropy. The input file contained 358 hadrons. 

The chemical potential depends on the strangeness, baryon number and charge of 
each particle, as well as the values of the corresponding potentials. Charge inde- 
pendence has been assumed (/xq = 0), which essentially disregards the last column 
of input. The strangeness potential is fixed by (jl.5J) . as the total strangeness of 
the system must be zero. 
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This leaves two free parameters: temperature T and baryonic chemical potential 
11b- By analysis of past heavy ion collisions, and applying the thermal model to the 
particle ratios observed, one can fit the parameters T, ji B and the other potentials 
to get the best agreement with experiment. These values correspond to the system 
parameters at chemical freeze-out. 

o.is | , , , , , , , , 




0.0-1 1 ■ 1 1 1 1 1 1 1 1 

0.1 0.2 0..* 0.-1 0.5 0.0 0." O.S 0." 
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Figure 3.1: Relationship between the temperature T and chemical potential fis at 
chemical freeze-out, as found by observing heavy ion collisions 



There are 25 data points shown in Fig 3.1, largely concentrated at T w 0.165GeV. 
The points furthest to the right correspond to the lowest energy collisions, and as 
the beam energies increase, the temperature T at freeze-out increases, while \xb 
decreases. As the baryonic chemical potential can never become negative (always 
excess baryons), and the temperature at the highest energy collisions appears to 
tend to an almost constant value, it is expected that when LHC becomes opera- 
tional, the parameters will have values fi B ~ + , and T ps 0.17 GeV. 

In all further calculations, the 25 data points describing the freeze-out curve are 
used. All plots shown are always against one of T or the other being implicit. 
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For all the points on the T-/i# curve, it has been pointed out by Cleymans and 
Redlich jl] that the average energy per hadron at the chemical freeze-out is almost 
constant throughout. In particular - w 1 GeV. 

Further, it has been shown by Cleymans 3j, that the baryonic potential can be 
fitted to the beam energy using a rather simple formula: 

~ 1 i a riu ( 3 - X ) 
1 + y/s/b 

with a fa 1.27 GeV, and b w 4.3 GeV. The temperature dependence on the beam 
energy can be implemented using - sa 1 GeV. Plots of the temperature and bary- 
onic potential as a function of \/s are provided in [3| here I only reproduce the 
curve for due to the unavailability of an analytic form of the T dependence. 
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Figure 3.2: Value of baryonic chemical potential \x& as a function of beam energy 
y/s, with parameter results from four experiments put it. 



Various features of the thermal model are calculated (next sections) in terms of T 
and fiB, and to find the corresponding beam energy \fs for a given (T, fin), the 
above relationship will be used. 
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3.1 Determining us 

The strangeness chemical potential [is is uniquely defined by T and [is, but there 
it however no analytic way of determining it directly. It needs to be solved numer- 
ically. To do so, it is essential to first know approximately what value is expected. 

3.1.1 First estimate 

Every strange particle has a corresponding anti-particle, which has exactly oppo- 
site strangeness, while identical mass and degeneracy factor. They have opposite 
charges, but with [iq = this is irrelevant. In the current discussion, a meson's 
chemical potential depends only on the strangeness of the particle, so of the two 
groups: strange mesons and anti-strange mesons, the one with the greater chemical 
potential will exist in large quantities. 

For baryons however, there is [Ib, which is always positive. Hence there will always 
be more baryons than anti-baryons, and in the baryon section of the fireball, there 
will be an excess of strange quarks. 

To keep total strangeness zero, there must be an excess of anti-strange quarks in 
the mesons. Hence /is-S(s) > 0. As strangeness of s is +1, we deduce that [is > 0. 

For an upper bound, I presumed that [is will scale not faster than linearly with 
[iB, and have taken [i s < \^b- Although I have no analytical way of showing that 
this will suffice, numerically it was shown to work. 

Combining two bounds, the original interval was taken as: < [is < \[Ib- 

3.1.2 Numerics 

The momentum distribution function was originally mentioned in the first chapter. 
It is defined as: 




-i -i 



fiip) 



e 



+ e 




Then the number of particles of species i that will be produced is: 




(3.3) 



Combining this with ()1.5|) . for a given T and [is' 




(3.4) 
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The problem has been transformed into one of finding a zero of the function T{[is)- 
Starting with the bounds < us < §/Ab, it turns out that the T is monotonically 
growing over this interval (^(0) < and jF(|/i B ) > 0), and an iterative solution 
(such as the bisection method) can be implemented to find a zero of T . 

Importantly here, each evaluation of T requires the evaluation of almost 200 in- 
tegrals (one only needs to find rii for strange particles) of a non-standard form. 
This can be done using the Simpson method, but that is extremely time consum- 
ing, and the program takes a very long time to find a value of us to a sufficient 
accuracy. I have calculated rii for a few typical parameter sets, and have found 
that for the 7r's the Simpson and Laguerre methods differ by around 20 ppm, while 
for all other particles the absolute difference was of order 10~ 8 or less. As pions 
are non-strange, all the necessary integrals can be computed using the very fast 
Laguerre method (15 points), and the bisection algorithm can be implemented. 

3.1.3 Results 
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Figure 3.3: Relationship between the strangeness us and baryonic potential ixb, 
calculated numerically along the freeze-out curve 
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The result is rather neat, although not too surprising. At very low when the 
temperature is almost constant, the strangeness potential goes linearly with 
As the baryonic potential increases, so does /is, but more slowly, and the plot 
curves downwards. The little twist at the last point is rather unexpected, possibly 
due to uncertainty in T. 

0.2-1 | , , , , , , , , 
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Figure 3.4: Ratio of strangeness to baryonic potential (i^s/^b) dependance on /is 



This plot is essentially the same as Fig 3.3, although it does show more clearly the 
linear dependence of //£ on \xb as the baryonic chemical potential is small, together 
with the fall-off for larger values of \xb (corresponding to lower temperatures). The 
limiting value for the ratio of the chemical potentials is: 

lim — « 0.236 (3.5) 

For a QGP, there are no hadrons, and for total strangeness to be zero, fi(s) = n(s), 
and as the strange quark has baryon number | and strangeness —1 (opposites signs 
for s), it is required that fj,s — §/Ab in a QGP. The limiting value found does not 
correspond to that for a QGP, as it is calculated for a hadron gas at freeze-out, 
which differs from the hadron gas at formation (when QGP disassociates), which in 
turn which differs from QGP through a phase transition, breaking limit continuity. 
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3.2 Entropy in a Hadron Gas 



Once the strangeness potential fi$ has been determined, the chemical potential for 
each hadron becomes known, and the integrals can be computed. Recalling that 
the entropy density for one hadron (|2.8jl can be written as: 

s * = / [-fi(x)lnf i (x)T(lTfi{x))HlTfi(x))]x 2 dx 

^ Jo 

it is expected that Sj will eventually go as T 3 for sufficiently large temperatures. 
Also, as the T 3 term is common for all hadrons at any given temperature, one can 
calculate the integrals without it, and insert it later. For this reason, I define the 
normalized entropy density for hadron i: 

h=^r [-fi(x) In fi(x) T (1 T fi(x)) ln(l =F fi(x))} x 2 dx 
^ Jo 

Then the entropy density of each group of hadrons can be written as: 



sm 


— T 3 Sm 


= t 3 £ h 






M esons 


SB 


= T*S B 








Baryons 


sr 


= T 3 S T 


= t 3 £ u 






Hadrons 



which defines the group normalized entropy densities <Sm, <Sb an d Sr- 



3.2.1 Meson-Baryon distribution 

The normalized entropy density fa was evaluated for each hadron. This was done 
using both the Simpson and Gauss-Laguerre integration techniques. It was ob- 
served, that although for several particles there was a "large" discrepancy (up to 
0.5%) between the two methods, these occurred only for the heaviest particles, 

where \j x 2 + ~ x is a very bad approximation for all but very large T. How- 
ever, the hadron gas is always dominated by particles of lowest mass. In the case of 
mesons these are the 7T, n, p, K's. For baryons it is p, n and light resonances, and 
for all of the mentioned particles, the integration methods do not differ by more 
that 10 ppm. If dealing with the (normalized) entropy density for mesons, baryons 
and the total, the two methods give very similar results (the relative differences 
never exceeded 20 ppm), which cannot be separated on a graph. 

The T-fiB relationship has been shown earlier (Fig 3.1), so all plots from here 
onwards will be against temperature. 
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It is expected that for systems with low energy (and temperature), the hadron gas 
will be dominated by the protons and neutrons (baryons), and that the mesons will 
exist only in small quantities. Here, I put anti-baryons together with baryons. For 
very high temperatures (as observed at RHIC) the system is however dominated 
by pions, and the entropy is expected to be mostly among the mesons. This is 
indeed observed: 

i-i , , , , , , , 
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Figure 3.5: Calculated normalized entropy density S for mesons, baryons, and 
total. The baryon and meson entropies are seen to cross when T ps 140MeV. 



The plot contains many interesting features. For low temperatures, as was ex- 
pected, the baryon entropy is far higher than that for mesons. Indeed, it appears 
to be far greater than the entropy at higher temperatures - this is however not true. 
The plot shows divided by the cube of the temperature. The baryon and meson 
curves intersect very near to T = 140MeV. This corresponds to /ib ~ 0.406 GeV, 
and y/s ~ 9.3 GeV. This falls within the domain of the Alternating Gradient Sys- 
tem (AGS) at the Brookhaven National Lab. For the high temperatures (lots of 
points, as T seems to approach a constant value as //# — > 0), the mesons dominate 
as expected, while the baryons continue to contribute substantially less for the 
increasing values of T. 
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The total normalized entropy density Sr exhibits a very neat behaviour. For 
T < lOOMeV, it decreases almost hyperbolically as the baryon contribution goes 
down. It then flattens out, has a minimum at T = 130MeV, and then increases 
very slowly. For the high T, the value of Sr is almost constant, as if it was 
approaching a limit. In that region Sr ~ 6.15 

This limiting behaviour is not totally surprising. It was expected that as T > m 
for the hadrons, sr ~ T 3 , which is the same as having Sr constant. However this 
limiting behaviour is observed already at T ~ 130MeV, which corresponds to the 
pion mass, and is a lot less than the masses of all the baryons. 

For completeness, one may also look at a plot of the proper entropy densities s 
and their dependance on temperature: 
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Figure 3.6: Calculated entropy density s for mesons, baryons, and the total. 



Here it can be clearly seen that the baryon entropy density, unlike Sq does not 
decrease initially. It increases (as one would expect), although slower that T 3 . It 
does however reach a maximum at T = 150MeV, when \xb starts decreasing fast 
relative to the increase in T. For mesons, the entropy density simply increases 
very fast with T. A power law fit gives sm ~ T 5 - 29±a04 . 
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3.2.2 Baryon - Anti-baryon distribution 



Although this was not not originally planned, it was an automatic consequence 
of doing the baryon-meson distributions, and though it does not give any new re- 
sults, it is a simple check on the intuitive expectations for a hadronic gas. At low 
temperatures, baryons are favoured greatly over anti-baryons. At higher temper- 
atures (and beam energies), hb — ► 0, and anti-particles are produced in sizeable 
quantities. Although there is still a bias towards baryons, it is far smaller. 

To do the plot, I originally planned to show Sq and Sq, but that plot is not clear. 
Recall that S for baryons plus anti-baryons is maximal for low T, and in that 
region S-g ~ 0. Only at high T will the anti-baryons contribute significantly, and 
so I plotted the proper entropy densities sq and Sg. 
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Figure 3.7: Calculated entropy densities sq and s^- 



Although maybe more pronounced at high T, the plot is as expected. The anti- 
baryons always contribute less than the baryons, but as hb ~^ 0, they do limit 
to the same value as the baryons. drops sharply as T 165.5 MeV remains 
almost constant while \ib decays almost linearly in the same region. 
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3.2.3 Entropy of strange particles 



The distribution of the entropy among the mesons, baryons and anti-baryons has 
been demonstrated. I have mentioned that for mesons, the entropy contributions 
come mainly from the n, i] and K particles, while baryon contributions come 
mainly from nucleons and light resonances. There are also lots of strange par- 
ticles with slightly higher masses, which for higher T are expected to contribute 
significantly. This does indeed happen: 
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Figure 3.8: Fractional contribution to the entropy of mesons/baryons/total from 
strange particles (the = ss is considered unflavoured) 



The result is rather interesting. For the mesons, even at low T, there is a little 
contribution from the kaons, and this rises steadily with increasing T. For the 
baryons, at low T there are almost no strange particles (as expected), but for higher 
T this rises very sharply, and for T > 160MeV, the strange particles (A, E etc.) 
dominate. This may be accidental, but the meson and baryon contributions cross 
around 0.5, corresponding to equal strange and non-strange contributions. 
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3.3 Entropy in a QGP 



The calculations of entropy in a QGP are simpler than those in a hadron gas, as 
there are far fewer "particles". In this case, the only contributors to the entropy 
are the quarks, anti-quarks and gluons. As for the hadrons, I only used unflavoured 
and strange quarks (u, d, s and the anti-quarks). Gluons can be treated in the 
same fashion. 

For baryon numbers, one takes B(q) = | and B(q) = — |. Charge independence 
gives jiQ = 0. To ensure that strangeness is conserved, /i(s) = /x(s), and hence 
Us — \hb- This deals with the three chemical potentials. 

The u and d quarks are assumed massless, and the strange quark was given a 
mass of 150MeV. Gluons are massless, and have no charge, baryon number or 
strangeness, so fj,(g) = 0. For spin degeneracies g s , I included colour, so g s (q) = 6 
and g s (g) = 8x2 = 16. This is all the information required to do the calculations. 
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Figure 3.9: Normalized entropy densities S for a QGP. The total is plotted, to- 
gether with the individual contributions from the quarks, anti-quarks and gluons 



25 



I had taken the T~hb set of parameters, and calculated the entropy density that 
would result from a QGP at those values (although jj, s is different). The result: 

As can be seen in Fig 3.9, the plot for the total entropy density is very similar to 
that of the total entropy density of a hadronic gas. Starting very large, it decreases 
almost hyperbolically, and then levels off at a value - presumably it will limit to a 
constant. There is however no little dip as in the HG case, and for large T, while 
seemingly limiting to a contant, the normalized density S is slowly decreasing for a 
QGP, whereas it was slighly increasing in the case of a HG. The expected limiting 
value for the QGP being S Qg -p m 20.45. 

Next, it is observed that the gluon normalized entropy density does not depend on 
temperature. This is exactly what is expected, as the gluons have no mass, and 
no chemical potential, so all the T factors drop out. For T > 165.5MeV, I found 
that S(g) > S(q) > S(q). It is required that the quarks dominate the anti-quarks, 
but at high temperatures the values become very close, and get exceeded by the 
gluon entropy. This was expected, although maybe only at higher T. 

The strangeness contribution to the total was also calculated: 
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Figure 3.10: Contribution of total entropy due to s and s quarks 
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The plot shown in Fig 3.10 does have a similar shape to the strangeness contri- 
bution of mesons in a hadron gas. There is however a large difference in that for 
mesons, the strange particle contribution went to a little more than a half, whereas 
here the strange (anti-strange) quark contributions are never nore than 0.21 of the 
total. This can be easily explained - for a HG, the strange (anti-strange) quarks 
can bond with an u or d quark to form light kaons. In QGP, strange quarks exist 
alone, and as they are heavy (m(s) ~ 150MeV), their entropy will be exceeded by 
the u quarks, the d quarks, AND the gluons at high energies. So the upper bound 
for s, s contribution is | = 0.25. Although this is not reached, the contribution is 
in the expected range. 



3.4 Comparison of HG and QGP 

The limiting value for the normalized entropy density in a QGP has been calculated 
to be a little more than three times greater than the corresponding value in a 
hadron gas. This is true approximately true for all T. A comparison plot for the 
total entropy densities s of the two phases is presented: 
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Figure 3.11: Total entropy densities s for a Hadron Gas and for a QGP 
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From Fig 3.11 one can clearly see that the entropy in a QGP is higher by a factor of 
between 3 and 5 than the entropy in a hadron gas with the same temperature and 
baryon chemical potential. This is due to the QGP consisting of essentially free 
quarks and gluons, allowing for far greater degrees of freedom, and hence directly 
increasing the entropy of the system. 

It is expected that the transition between a hadron gas and a QGP will be through 
a phase transition. If one could take a hadron gas and increase the temperature 
until a QGP formed, the entropy would originally follow the HG curve, and at 
some critical temperature "jump" to the QGP curve. As the entropy curves differ 
by a an almost constant factor (~ 3), a jump discontinuity is expected, signifying 
a first-order phase transition. 

3.5 Conclusion 

For an ideal hadron gas model of heavy ion collisions, the entropy contributions 
have been calculated at freeze-out for a set of parameter values satisfying the 
Cleymans-Redlich freeze-out criterion. The entropy densities have been calculated 
for mesons and baryons as well as the total (Fig 3.5 and Fig 3.6). It has been 
demonstrated that at low temperatures (energies) the entropy of the system is 
baryon dominated, while at high T the mesons dominate. The baryon-meson 
cross-over occurs for T ~ 140 MeV, which corresponds to ^/s pa 9.3 GeV. This 
falls within the domain of the AGS at Brookhaven. 

For each of the considered parameter values (T and hb) the entropy density of a 
corresponding QGP has been calculated and presented (Fig 3.11). Due to "free" 
quarks in a QGP and extra degrees of freedom, the entropy in each corresponding 
QGP is significantly higher than that of the HG. This shows that for a transi- 
tion from a hadron gas to a QGP, the entropy would have a jump discontinuity, 
characteristic of a first order phase transition. 

At very low temperatures T pa 0.05 GeV , the strange particle contributions to- 
wards the entropy (in both HG and QGP) are almost negligible. As the energies 
increase, the strange particle contributions increase monotonically with T, up to 
55% in a hadron gas, and almost 25% in a QGP. 

Also, for a hadron gas the baryon - anti-baryon dependence was investigated. 
At medium and low temperatures (T < 140 MeV), the contributions from anti- 
baryons remain negligible. Only at very high energies (such as at RHIC), when 
T pa 165 MeV, the multiplicities of anti-baryons become significant, and the ratio 
2 pa 0.8 for the highest T data points. For LHC it is expected that 2 — > \~ . 
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